
#####################################################
#
#               Code to replicate 
#  Affluence and influence in a social democracy 
#   
#    Ruben Berge Mathisen (University of Bergen)
#
#       Tables and Figures for Appendix
#                February 2022
#
#   Original analyses conducted in R version 3.5.1 
#
#####################################################



#####################################################
######## Packages and data prep #####################
#####################################################

library(car)
library(dplyr)
library(tidyr)
library(readxl)
library(kableExtra)
library(ggplot2)
library(ggrepel)
library(stringr)
library(stargazer)
library(cowplot)
library(gridExtra)

d <- read_excel("Data for Affluence and influence in a social democracy Feb 2022.xlsx")

# Logit-transform variables
logitTransform <- function(p) { log(p/(1-p)) }

rec <- function(x) {
  x <- ifelse(x>=1,0.9999,ifelse(x<=0,0.0001,x))
  logitTransform(x)  
}
d <- d %>% mutate_at(vars(pred_i90e90:pred_i10e10),funs(rec))
d$all_logit <- logitTransform(d$all)
d$inc_p90_logit <- logitTransform(d$inc_p90)
d$inc_p70_logit <- logitTransform(d$inc_p70)
d$inc_p50_logit <- logitTransform(d$inc_p50)
d$inc_p30_logit <- logitTransform(d$inc_p30)
d$inc_p10_logit <- logitTransform(d$inc_p10)
d$edu_p90_logit <- logitTransform(d$edu_p90)

d$men_logit <- logitTransform(d$men)
d$women_logit <- logitTransform(d$women)

d_all <- d
# Remove half-adopted policies 
d <- d %>% filter(outcome!=0.5)
# REMOVE POLICIES THAT REQUIRE CONSTITUIONAL CHANGE
d <- d %>% filter(constitutional_change!="yes")
# Remove academic
d_aca <- d 
d <- d %>% filter(academic_survey==0)


#####################################################
######## Functions ##################################
#####################################################

# Functions for calulcating table and figure quantities
ests <- NULL
eff<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  ests <- data.frame(rbind(ests,round(summary(m)$coefficients[2,1],2)))
  p <- summary(m)$coefficients[2,4]
  se <- round(summary(m)$coefficients[2,2],2)
  p1 <- car::recode(p,"lo:0.009999='***'; 0.01:0.049999='**'; 0.05:0.099999='*'; 0.1:hi=''")
  result <- paste0(ests, p1)
} 
eff_se<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  ests <- data.frame(rbind(ests,round(summary(m)$coefficients[2,1],2)))
  se <- round(summary(m)$coefficients[2,2],2)
  result <- paste0(ests," (",se,")")
} 
eff_num<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  ests <- data.frame(rbind(ests,round(summary(m)$coefficients[2,1],2)))
  ests
} 
se<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  se <- round(summary(m)$coefficients[2,2],2)
  se
} 
p<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  p <- summary(m)$coefficients[2,4]
  substring(as.character(ifelse(p<0.001,"<<.001", round(p,3))),2)
} 
sup20<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  round(predict(m,newdata=data.frame(x=logitTransform(0.2)),type="response"),2)
} 
sup80<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  round(predict(m,newdata=data.frame(x=logitTransform(0.8)),type="response"),2)
} 
fact<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  diff<- predict(m,newdata=data.frame(x=logitTransform(0.8)),type="response")/predict(m,newdata=data.frame(x=logitTransform(0.2)),type="response")
  round(diff,1)
} 
n<-function(x){
  m<-glm(outcome~x,data=d1,family="binomial")
  length(m$residuals)
} 


#####################################################
##### Tables and Figures (same order as in paper) ###
#####################################################

# Table A1 
rec <- function(x) round(x,2)*100
d1 <- d %>% 
  group_by(topic) %>% 
  arrange(-diff_inc90_inc10) %>%
  select(year,policy,inc_p10,inc_p50,inc_p90,inc90_minus_inc10,outcome) %>%
  mutate_at(vars(inc_p10:inc90_minus_inc10),funs(rec)) %>%
  mutate(inc90_minus_inc10 = ifelse(inc90_minus_inc10>0,paste0("+",inc90_minus_inc10),inc90_minus_inc10),
         outcome=ifelse(outcome==1,"YES","NO")) %>%
  slice(1:10) 

kable(d1[,-1],booktabs=T,align="l", caption = "Most constested policy proposals across income levels, by policy area",
      col.names = c(" ","Policy","Poor (P10)","Middle (P50)","Affluent (P90)","P90-P10", "Adopted after 4 years?")) %>%
  pack_rows("Economy",1,10) %>%
  pack_rows("Moral",11,20) %>%
  pack_rows("Foreign/security",21,30) %>%
  pack_rows("Other",31,40) %>% 
  column_spec(1,width = "0.5in") %>%
  column_spec(3:7,width = "0.4in") %>%
  add_header_above(c(" "=2, "Support among 
                     income groups (%)" = 3, " "=2)) %>%
  kable_styling(font_size = 9) %>%
  kable_styling(latex_options = "HOLD_position")

# Table A2 
d1 <- d_all %>% group_by(topic,topic_specific) %>% tally() %>% arrange(topic,-n)
d2 <- d %>% group_by(topic,topic_specific) %>% tally() %>% arrange(topic,-n)
d3 <- merge(d1,d2,by=c("topic","topic_specific"),all.x = T)
d3$n.y <- ifelse(is.na(d3$n.y),0,d3$n.y) 
d3$var <- paste0(d3$topic_specific," (",d3$n.x,"/",d3$n.y ,")")
d3 <- d3 %>% dplyr::select(topic,var)
d3_eco <- d3 %>% filter(topic=="Economy") %>% dplyr::select(var) %>% dplyr::mutate(row=as.numeric(row.names(.)))
d3_mor <- d3 %>% filter(topic=="Moral") %>% dplyr::select(var) %>% dplyr::mutate(row=as.numeric(row.names(.)))
d3_for <- d3 %>% filter(topic=="Foreign/Security") %>% dplyr::select(var) %>% dplyr::mutate(row=as.numeric(row.names(.)))
d3_oth <- d3 %>% filter(topic=="Other") %>% dplyr::select(var) %>% dplyr::mutate(row=as.numeric(row.names(.)))
d4 <- merge(d3_eco,d3_mor,by="row", all.y=T)
d4 <- merge(d4,d3_for,by="row", all.x=T)
d4 <- merge(d4,d3_oth,by="row", all.x=T)
kable(d4[,-1],
      booktabs=T,
      align="l",
      row.names = F,
      linesep = "",
      col.names = c("Economic policy (162/110)","Moral policy; value issues (212/164)","Foreign/Security policy (134/68)","Other issues (95/55)"),
      caption="Coding of policy areas (n in full dataset / n in dataset used in letter)") %>% 
  kable_styling(font_size = 10) %>%
  kable_styling(latex_options = "HOLD_position") %>%
  column_spec(1,width = "1.8in") %>%
  column_spec(2,width = "1.5in") %>%
  column_spec(3,width = "1.7in") %>%
  column_spec(4,width = "1.5in")


# Table A3
d1 <- d %>% filter(diff_inc90_inc10>0.08&topic=="Economy")
tab1<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90a=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
d1 <- d %>% filter(diff_inc90_inc10>0.08&topic=="Moral")
tab2<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90b=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
tab <- cbind(tab1,tab2)
d1 <- d %>% filter(diff_inc90_inc10>0.04&topic=="Foreign/Security")
tab3<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90b=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
d1 <- d %>% filter(diff_inc90_inc10>0.04&topic=="Other")
tab4<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90b=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
tab <- cbind(tab1,tab2,tab3,tab4)
tab1 <- data.frame(t(tab))
row.names(tab1) <- c("Poor (P10)",
                     "Affluent (P90) ",
                     "Poor (P10) ",
                     "Affluent (P90)  ",
                     "Poor (P10)  ",
                     "Affluent (P90)   ",
                     "Poor (P10)   ",
                     "Affluent (P90)    ")
kable(tab1[-3],booktabs=T, align="c",linesep = '\\addlinespace', col.names = c("Effect (logit coefficient)","Standard error", "Predicted probability of policy change at 20% support","Predicted probability of policy change at 80% support","Relative change in probability","N"),
      caption = "Policy reponsiveness when preferences diverege, by policy area.") %>%
  column_spec(1,width="1.4in") %>%
  column_spec(2,width="0.5in") %>% 
  column_spec(3,width="0.55in") %>% 
  column_spec(4:5,width="0.8in") %>% 
  column_spec(6,width="0.5in") %>% 
  column_spec(7,width="0.4in") %>% 
  pack_rows("Economy",1,2) %>%
  pack_rows("Moral",3,4) %>%
  pack_rows("Foreign/Security",5,6) %>%
  pack_rows("Other",7,8) %>%
  footnote("Bivariate logistic regression models (rows). The dependent variable is a dichotomous measure of whether or not the policy change was adopted within four years of the time of the survey question. Included are policy proposals where preferences diverege between the 90th and 10th income percentiles. To ensure an acceptable sample size in each domain (around 50 or more where possible), the cutoff-point for inclusion was set to 8 points for economic and moral issues, and 4 points for foreign policy and 'other' issues (due to the lower n on the latter two topics). The lower threshold for the latter two could pose a challenge for comparison if it showed equal responsiveness, however, since it already shows highly unequal responsiveness, we can be fairly certain that this would not change much with a higher threshold. The same, of course, applies to moral issues. The effect size for the poor on economic policy remained essentially the same at higher cut-off points (b=0.62, se=0.33 at 10p; b=0.58, se=0.37 at 12p). The results are also robust to an alternative model specification, as reported in Table A4 below.*p<0.1; **p<0.05; ***p<0.01",threeparttable = T) %>%
  kable_styling(latex_options = "HOLD_position")

# Table A4
m1<-summary(lm(outcome~inc90_minus_inc10+all,data=d%>%filter(topic=="Economy")))
m2<-summary(lm(outcome~inc90_minus_inc10+all,data=d%>%filter(topic=="Moral")))
m3<-summary(lm(outcome~inc90_minus_inc10+all,data=d%>%filter(topic=="Foreign/Security")))
m4<-summary(lm(outcome~inc90_minus_inc10+all,data=d%>%filter(topic=="Other")))
tab <- data.frame(area=c("Affluent minus poor",
                         "Affluent minus poor",
                         "Affluent minus poor",
                         "Affluent minus poor"),
                  coef=round(c(m1$coeff[2,1],
                               m2$coeff[2,1],
                               m3$coeff[2,1],
                               m4$coeff[2,1]),2),
                  se=round(c(m1$coeff[2,2],
                             m2$coeff[2,2],
                             m3$coeff[2,2],
                             m4$coeff[2,2]),2),
                  p=round(c(m1$coeff[2,4],
                            m2$coeff[2,4],
                            m3$coeff[2,4],
                            m4$coeff[2,4]),3),
                  n=c(length(m1$residuals),
                      length(m2$residuals),
                      length(m3$residuals),
                      length(m4$residuals)),
                  stars=c(" ",
                          "***",
                          "***",
                          "***"))
tab$coef <- paste0(tab$coef,tab$stars)
kable(tab[,-6],booktabs=T, align="c",linesep = '\\addlinespace',
      col.names = c(" ", 
                    "Effect (OLS coefficient)",
                    "Standard error",
                    "p-value",
                    "N"),
      caption="Alternative specification for estimating unequal responsiveness by policy area") %>%
  pack_rows("Economy",1,1) %>%
  pack_rows("Moral",2,2) %>%
  pack_rows("Foreign/Security",3,3) %>%
  pack_rows("Other",4,4) %>%
  column_spec(1,width="1.7in") %>%
  column_spec(2:4,width="0.8in") %>%
  footnote("OLS regression models (rows). Alternative specification, as proposed by Schakel and Hakhverdian (2020, 154-155). The dependent variable is a dichotomous measure of whether or not the policy change was adopted within four years of the time of the survey question. The independent variable 'Affluent minus poor' was calculated by taking (% support of P90) minus (% support of P10). A positive coefficient means that policy is biased towards the preferences of the affluent, while a negative coefficient means that policy is biased towards the preferences of the poor. All models include controls for overall policy support. See Schakel and Hakhverdian (2020) for details about this method. *p<0.1; **p<0.05; ***p<0.01",threeparttable = T) %>%  
  kable_styling(latex_options = "HOLD_position")


# Figure B1
## 90 vs 10
d2 <- d
sequence <- seq(0,0.15,0.01)
est <- NULL
pval <- NULL
obs <- NULL
for(i in sequence){
  d1 <- d2 %>% filter(diff_inc90_inc10>i)
  fit = summary(glm(outcome~inc_p90_logit,data=d1,family="binomial"))
  obs <- rbind(obs,nrow(subset(d1,!is.na(outcome))))
  est <- data.frame(rbind(est,fit$coefficients[,1]))
  est$percentile <- "90th"
  pval <- data.frame(rbind(pval,fit$coefficients[,4]))
  cut90 <- cbind(est,pval,obs)[,-c(1,4)]
  names(cut90)[1] <- "coef"
  names(cut90)[3] <- "pval"
}
sequence <- seq(0,0.15,0.01)
est <- NULL
pval <- NULL
obs <- NULL
for(i in sequence){
  d1 <- d2 %>% filter(diff_inc90_inc10>i)
  fit = summary(glm(outcome~inc_p10_logit,data=d1,family="binomial"))
  obs <- rbind(obs,nrow(subset(d1,!is.na(outcome))))
  est <- data.frame(rbind(est,fit$coefficients[,1]))
  est$percentile <- "10th"
  pval <- data.frame(rbind(pval,fit$coefficients[,4]))
  cut10 <- cbind(est,pval,obs)[,-c(1,4)]
  names(cut10)[1] <- "coef"
  names(cut10)[3] <- "pval"
}
cut10$cutoff <- sequence
cut90$cutoff <- sequence
cuts <- rbind(cut90,cut10)
cuts$sig <- car::recode(cuts$pval,"lo:0.01='***'; 0.01:0.05='**'; 0.05:0.1='*'; 0.1:hi='ns'")
cuts$hjust <- ifelse(cuts$percentile=="90th", 0, 1)
cuts$cutoff_n <- paste0(cuts$cutoff,"
                        (",cuts$obs,")")
cuts9010 <- cuts
cuts9010$comp <- "Alternative thresholds for preference divergence (n): Poor v Affluent"
## 90 vs 50
d2 <- d
sequence <- seq(0,0.13,0.01)
est <- NULL
pval <- NULL
obs <- NULL
for(i in sequence){
  d1 <- d2 %>% filter(diff_inc90_inc50>i)
  fit = summary(glm(outcome~inc_p90_logit,data=d1,family="binomial"))
  obs <- rbind(obs,nrow(subset(d1,!is.na(outcome))))
  est <- data.frame(rbind(est,fit$coefficients[,1]))
  est$percentile <- "90th"
  pval <- data.frame(rbind(pval,fit$coefficients[,4]))
  cut90 <- cbind(est,pval,obs)[,-c(1,4)]
  names(cut90)[1] <- "coef"
  names(cut90)[3] <- "pval"
}
sequence <- seq(0,0.13,0.01)
est <- NULL
pval <- NULL
obs<- NULL
for(i in sequence){
  d1 <- d2 %>% filter(diff_inc90_inc50>i)
  fit = summary(glm(outcome~inc_p50_logit,data=d1,family="binomial"))
  obs <- rbind(obs,nrow(subset(d1,!is.na(outcome))))
  est <- data.frame(rbind(est,fit$coefficients[,1]))
  est$percentile <- "50th"
  pval <- data.frame(rbind(pval,fit$coefficients[,4]))
  cut50 <- cbind(est,pval,obs)[,-c(1,4)]
  names(cut50)[1] <- "coef"
  names(cut50)[3] <- "pval"
}
cut50$cutoff <- sequence
cut90$cutoff <- sequence
cuts <- rbind(cut90,cut50)
cuts$sig <- car::recode(cuts$pval,"lo:0.01='***'; 0.01:0.05='**'; 0.05:0.1='*'; 0.1:hi='ns'")
cuts$hjust <- ifelse(cuts$percentile=="90th", 0, 1)
cuts$cutoff_n <- paste0(cuts$cutoff,"
                        (",cuts$obs,")")
cuts9050 <- cuts
cuts9050$comp <- "Alternative thresholds for preference divergence (n): Middle v Affluent"
## Combine
cuts <- rbind(cuts9010,cuts9050)
cuts$comp <- factor(cuts$comp,levels(factor(cuts$comp))[c(1,2)])
cuts$label <- ifelse(cuts$cutoff==0.02&cuts$percentile=="90th","Affluent",
                     ifelse(cuts$cutoff==0.03&cuts$percentile=="10th","Poor",
                            ifelse(cuts$cutoff==0.03&cuts$percentile=="50th","Middle",NA)))
ggplot(subset(cuts,!is.na(cutoff_n)),aes(x=cutoff_n,
                                         y=coef,
                                         group=percentile,
                                         color=percentile,
                                         shape=percentile)) +
  geom_line() +
  geom_point() +
  labs(x=NULL,y="Opinion-policy link (logistic coefficient)\n", 
       group="Income percentile", 
       linetype="Income percentile") +
  theme_bw() +
  scale_y_continuous(limits=c(-0.4,1.2)) +
  geom_hline(yintercept=0,color="grey60") +
  theme(legend.position = "none",
        panel.grid = element_blank(),
        text = element_text(size = 15)) +
  geom_text(aes(label=sig,hjust=hjust), size=5, vjust=0) +
  facet_wrap(~comp,ncol=1,scales="free_x") +
  geom_text(aes(label=label),vjust=ifelse(cuts$percentile=="90th",-3.5,3.5))


# Table C1
d1 <- d %>% filter(diff_inc90_inc10>0.1&year<1985)
tab1<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90a=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
d1 <- d %>% filter(diff_inc90_inc10>0.1&year>1984)
tab2<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90b=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
tab <- cbind(tab1,tab2)

tab1 <- data.frame(t(tab))
row.names(tab1) <- c("Poor (P10)",
                     "Affluent (P90) ",
                     "Poor (P10) ",
                     "Affluent (P90)")

kable(tab1[-3],booktabs=T, align="c",linesep = '\\addlinespace', col.names = c("Effect (logit coefficient)","Standard error", "Predicted probability of policy change at 20% support","Predicted probability of policy change at 80% support","Relative change in probability","N"),
      caption = "Policy reponsiveness when preferences diverege (>10p), by time period.") %>%
  column_spec(1,width="1.4in") %>%
  column_spec(2,width="0.5in") %>% 
  column_spec(3,width="0.55in") %>% 
  column_spec(4:5,width="0.8in") %>% 
  column_spec(6,width="0.5in") %>% 
  column_spec(7,width="0.4in") %>% 
  pack_rows("1966-1984",1,2) %>%
  pack_rows("1985-2014",3,4) %>%
  footnote("Bivariate logistic regression models (rows). The dependent variable is a dichotomous measure of whether or not the policy change was adopted within four years of the time of the survey question. *p<0.1; **p<0.05; ***p<0.01",threeparttable = T) %>%
  kable_styling(latex_options = "HOLD_position")


# Table C2
d1 <- d %>% filter(diff_inc90_inc10>0.1&share_dk>0.13)
tab1<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90a=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
d1 <- d %>% filter(diff_inc90_inc10>0.1&share_dk<0.13)
tab2<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90b=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
tab <- cbind(tab1,tab2)
tab1 <- data.frame(t(tab))
row.names(tab1) <- c("Poor (P10)",
                     "Affluent (P90) ",
                     "Poor (P10) ",
                     "Affluent (P90)")
kable(tab1[-3],booktabs=T, align="c",linesep = '\\addlinespace', col.names = c("Effect (logit coefficient)","Standard error", "Predicted probability of policy change at 20% support","Predicted probability of policy change at 80% support","Relative change in probability","N"),
      caption = "Policy reponsiveness when preferences diverege (>10p), by issue salience.") %>%
  column_spec(1,width="1.4in") %>%
  column_spec(2,width="0.5in") %>% 
  column_spec(3,width="0.55in") %>% 
  column_spec(4:5,width="0.8in") %>% 
  column_spec(6,width="0.5in") %>% 
  column_spec(7,width="0.4in") %>% 
  pack_rows("Low salience (% 'Don't knows' > 13)",1,2) %>%
  pack_rows("High salience (% 'Don't knows' < 13)",3,4) %>%
  footnote("Bivariate logistic regression models (rows). The dependent variable is a dichotomous measure of whether or not the policy change was adopted within four years of the time of the survey question. *p<0.1; **p<0.05; ***p<0.01",threeparttable = T) %>%
  kable_styling(latex_options = "HOLD_position")



# Table D1
d1 <- d_all %>% group_by(academic_survey) %>% dplyr::summarise(n = length(policy),
                                                               p10 = mean(inc_p10),
                                                               p50 = mean(inc_p50),
                                                               p90 = mean(inc_p90),
                                                               diff=mean(diff_inc90_inc10),
                                                               out=mean(outcome),
                                                               grad=mean(grad),
                                                               uni = 1-length(unique(policy))/length(policy)) %>%
  mutate_at(vars(p10:uni),funs(round(.,2))) 

d1 <- data.frame(t(d1))
d1[2,] <- as.character(d1[2,])
row.names(d1) <- c(" ","N (policy propals)", "Mean support P10","Mean support P50","Mean support P90","Mean absolute support distance (P10,P90)","Share of policies adopted","Share of gradual questions ('reduce', 'increase', 'expand', etc.)","Share of proposals asked about more than once")
kable(d1[-1,],booktabs=T,linesep = '\\addlinespace',col.names = c("Commercial survey","Academic survey"),align="c",caption="Differences and similarities between commercial and academic survey data")  %>%
  column_spec(1,width="2in") %>%
  column_spec(2:3,width="0.7in") %>%
  add_header_above(c(" "=1,"Data source"=2)) %>%
  kable_styling(latex_options = "HOLD_position")
  

# Table D2
d1 <- d_aca
tab<-data.frame(all=t(d1 %>% summarise_at(vars(all_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                inc30=t(d1 %>% summarise_at(vars(inc_p30_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                inc50=t(d1 %>% summarise_at(vars(inc_p50_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                inc70=t(d1 %>% summarise_at(vars(inc_p70_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                inc90=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
tab1 <- data.frame(t(tab))
row.names(tab1) <- c("All","P10","P30","P50","P70","P90")
kable(tab1[-3], align="c",booktabs=T,linesep = '\\addlinespace', 
      col.names = c("Effect (logit coefficient)",
                    "Standard error",
                    "Predicted probability of policy change at 20% support",
                    "Predicted probability of policy change at 80% support",
                    "Relative change in probability",
                    "N"),
      caption = "Policy responsiveness by income, when including the data from academic surveys (Norwegian Election studies 1969-2013, and the Norwegian Citizen Panel 2013-2014; not included in the data presented in the main maunscript).") %>%
  column_spec(1,width="0.9in") %>%
  column_spec(2:3,width="0.6in") %>% 
  column_spec(4:5,width="0.8in") %>% 
  column_spec(6,width="0.5in") %>% 
  column_spec(7,width="0.4in") %>% 
  pack_rows("Income percentile",2,6) %>%
  footnote("Bivariate logistic regression models (rows). The dependent variable is a dichotomous measure of whether or not the policy change was adopted within four years of the time of the survey question. *p<0.1; **p<0.05; ***p<0.01",threeparttable = T) %>%
  kable_styling(latex_options = "HOLD_position")


# Table D3
d1 <- d_aca %>% filter(diff_inc90_inc50>0.1)
tab1<-data.frame(inc50=t(d1 %>% summarise_at(vars(inc_p50_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90a=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
d1 <- d_aca %>% filter(diff_inc90_inc10>0.1)
tab2<-data.frame(inc10=t(d1 %>% summarise_at(vars(inc_p10_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90b=t(d1 %>% summarise_at(vars(inc_p90_logit),funs(eff,se,p,sup20,sup80,fact,n))))
tab <- cbind(tab2,tab1)
tab1 <- data.frame(t(tab))
row.names(tab1) <- c("Poor (P10)",
                     "Affluent (P90) ",
                     "Middle (P50)",
                     "Affluent (P90)")
kable(tab1[-3],booktabs=T, align="c",linesep = '\\addlinespace', col.names = c("Effect (logit coefficient)","Standard error", "Predicted probability of policy change at 20% support","Predicted probability of policy change at 80% support","Relative change in probability","N"),
      caption = "Policy responsiveness when preference diverege by more than 10 points, when including data from academic surveys.") %>%
  column_spec(1,width="1.4in") %>%
  column_spec(2,width="0.5in") %>% 
  column_spec(3,width="0.55in") %>% 
  column_spec(4:5,width="0.8in") %>% 
  column_spec(6,width="0.5in") %>% 
  column_spec(7,width="0.4in") %>% 
  pack_rows("Poor vs. affluent",1,2) %>%
  pack_rows("Middle vs. affluent",3,4) %>%
  footnote("Bivariate logistic regression models (rows). The dependent variable is a dichotomous measure of whether or not the policy change was adopted within four years of the time of the survey question. *p<0.1; **p<0.05; ***p<0.01",threeparttable = T) %>%
  kable_styling(latex_options = "HOLD_position")


# Figure D1
## 90 vs 10
d2 <- d_aca
sequence <- seq(0,0.15,0.01)
est <- NULL
pval <- NULL
obs <- NULL
for(i in sequence){
  d1 <- d2 %>% filter(diff_inc90_inc10>i)
  fit = summary(glm(outcome~inc_p90_logit,data=d1,family="binomial"))
  obs <- rbind(obs,nrow(subset(d1,!is.na(outcome))))
  est <- data.frame(rbind(est,fit$coefficients[,1]))
  est$percentile <- "90th"
  pval <- data.frame(rbind(pval,fit$coefficients[,4]))
  cut90 <- cbind(est,pval,obs)[,-c(1,4)]
  names(cut90)[1] <- "coef"
  names(cut90)[3] <- "pval"
}
sequence <- seq(0,0.15,0.01)
est <- NULL
pval <- NULL
obs <- NULL
for(i in sequence){
  d1 <- d2 %>% filter(diff_inc90_inc10>i)
  fit = summary(glm(outcome~inc_p10_logit,data=d1,family="binomial"))
  obs <- rbind(obs,nrow(subset(d1,!is.na(outcome))))
  est <- data.frame(rbind(est,fit$coefficients[,1]))
  est$percentile <- "10th"
  pval <- data.frame(rbind(pval,fit$coefficients[,4]))
  cut10 <- cbind(est,pval,obs)[,-c(1,4)]
  names(cut10)[1] <- "coef"
  names(cut10)[3] <- "pval"
}
cut10$cutoff <- sequence
cut90$cutoff <- sequence
cuts <- rbind(cut90,cut10)
cuts$sig <- car::recode(cuts$pval,"lo:0.01='***'; 0.01:0.05='**'; 0.05:0.1='*'; 0.1:hi='ns'")
cuts$hjust <- ifelse(cuts$percentile=="90th", 0, 1)
cuts$cutoff_n <- paste0(cuts$cutoff,"
                        (",cuts$obs,")")
cuts9010 <- cuts
cuts9010$comp <- "Alternative thresholds for preference divergence (n): Poor v Affluent"
## 90 vs 50
d2 <- d_aca
sequence <- seq(0,0.15,0.01)
est <- NULL
pval <- NULL
obs <- NULL
for(i in sequence){
  d1 <- d2 %>% filter(diff_inc90_inc50>i)
  fit = summary(glm(outcome~inc_p90_logit,data=d1,family="binomial"))
  obs <- rbind(obs,nrow(subset(d1,!is.na(outcome))))
  est <- data.frame(rbind(est,fit$coefficients[,1]))
  est$percentile <- "90th"
  pval <- data.frame(rbind(pval,fit$coefficients[,4]))
  cut90 <- cbind(est,pval,obs)[,-c(1,4)]
  names(cut90)[1] <- "coef"
  names(cut90)[3] <- "pval"
}
sequence <- seq(0,0.15,0.01)
est <- NULL
pval <- NULL
obs<- NULL
for(i in sequence){
  d1 <- d2 %>% filter(diff_inc90_inc50>i)
  fit = summary(glm(outcome~inc_p50_logit,data=d1,family="binomial"))
  obs <- rbind(obs,nrow(subset(d1,!is.na(outcome))))
  est <- data.frame(rbind(est,fit$coefficients[,1]))
  est$percentile <- "50th"
  pval <- data.frame(rbind(pval,fit$coefficients[,4]))
  cut50 <- cbind(est,pval,obs)[,-c(1,4)]
  names(cut50)[1] <- "coef"
  names(cut50)[3] <- "pval"
}
cut50$cutoff <- sequence
cut90$cutoff <- sequence
cuts <- rbind(cut90,cut50)
cuts$sig <- car::recode(cuts$pval,"lo:0.01='***'; 0.01:0.05='**'; 0.05:0.1='*'; 0.1:hi='ns'")
cuts$hjust <- ifelse(cuts$percentile=="90th", 0, 1)
cuts$cutoff_n <- paste0(cuts$cutoff,"
                        (",cuts$obs,")")
cuts9050 <- cuts
cuts9050$comp <- "Alternative thresholds for preference divergence (n): Middle v Affluent"
## Combine
cuts <- rbind(cuts9010,cuts9050)
cuts$comp <- factor(cuts$comp,levels(factor(cuts$comp))[c(1,2)])
cuts$label <- ifelse(cuts$cutoff==0.02&cuts$percentile=="90th","Affluent",
                     ifelse(cuts$cutoff==0.03&cuts$percentile=="10th","Poor",
                            ifelse(cuts$cutoff==0.03&cuts$percentile=="50th","Middle",NA)))
ggplot(subset(cuts,!is.na(cutoff_n)),aes(x=cutoff_n,
                                         y=coef,
                                         group=percentile,
                                         color=percentile,
                                         shape=percentile)) +
  geom_line() +
  geom_point() +
  labs(x=NULL,y="Opinion-policy link (logistic coefficient)\n", 
       group="Income percentile", 
       linetype="Income percentile") +
  theme_bw() +
  scale_y_continuous(limits=c(-0.2,1.9)) +
  geom_hline(yintercept=0,color="grey60") +
  theme(legend.position = "none",
        panel.grid = element_blank(),
        text = element_text(size = 15)) +
  geom_text(aes(label=sig,hjust=hjust), size=5, vjust=0) +
  facet_wrap(~comp,ncol=1,scales="free_x") +
  geom_text(aes(label=label),vjust=ifelse(cuts$percentile=="90th",-3.5,3.5))


# Table E1
d1 <- d %>% filter(diff_inc90_inc10>0.1 | diff_edu90_edu10>0.1) 
tab1<-data.frame(inc10a=t(d1 %>% summarise_at(vars(pred_i10e10),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc10b=t(d1 %>% summarise_at(vars(pred_i50e10),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc10c=t(d1 %>% summarise_at(vars(pred_i90e10),funs(eff,se,p,sup20,sup80,fact,n))))
d1 <- d %>% filter((diff_inc90_inc10>0.1 | diff_edu90_edu10>0.1)) 
tab2<-data.frame(inc50a=t(d1 %>% summarise_at(vars(pred_i10e50),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc50b=t(d1 %>% summarise_at(vars(pred_i50e50),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc50c=t(d1 %>% summarise_at(vars(pred_i90e50),funs(eff,se,p,sup20,sup80,fact,n))))
d1 <- d %>% filter((diff_inc90_inc10>0.1 | diff_edu90_edu10>0.1)) 
tab3<-data.frame(inc90a=t(d1 %>% summarise_at(vars(pred_i10e90),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90b=t(d1 %>% summarise_at(vars(pred_i50e90),funs(eff,se,p,sup20,sup80,fact,n))),
                 inc90c=t(d1 %>% summarise_at(vars(pred_i90e90),funs(eff,se,p,sup20,sup80,fact,n))))
tab <- cbind(tab1,tab2,tab3)
tab1 <- data.frame(t(tab))
row.names(tab1) <- c("Income P10",
                     "Income P50",
                     "Income P90",
                     "Income P10 ",
                     "Income P50 ",
                     "Income P90 ",
                     "Income P10  ",
                     "Income P50  ",
                     "Income P90  ")
kable(tab1[-3],booktabs=T, align="c",linesep = '\\addlinespace', col.names = c("Effect (logit coefficient)","Standard error", "Predicted probability of policy change at 20% support","Predicted probability of policy change at 80% support","Relative change in probability","N"),
      caption = "Policy responsiveness to income/education percentile.") %>%
  column_spec(1,width="1.4in") %>%
  column_spec(2,width="0.5in") %>% 
  column_spec(3,width="0.55in") %>% 
  column_spec(4:5,width="0.8in") %>% 
  column_spec(6,width="0.5in") %>% 
  column_spec(7,width="0.4in") %>% 
  pack_rows("Education P10",1,3) %>%
  pack_rows("Education P50",4,6) %>%
  pack_rows("Education P90",7,9) %>%  
  footnote("Bivariate logistic regression models (rows) for 9 combinations of education and income percentile. The support among an income/education combination is the logit-transformed imputed percentage of respondents favoring the policy change in that income/education combination. Included are policy proposals with minimum 10 point preference divergence between the 90th and 10th income percentiles or 90th and 10th education percentiles. The dependent variable is a dichotomous measure of whether or not the policy change was adopted within four years of the time of the survey question. *p<0.1; **p<0.05; ***p<0.01",threeparttable = T) %>%
  kable_styling(latex_options = "HOLD_position")


# Table G1
d1 <- d 
tab1<-data.frame(men=t(d1 %>% summarise_at(vars(men_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 women=t(d1 %>% summarise_at(vars(women_logit),funs(eff,se,p,sup20,sup80,fact,n))))
d1 <- d %>% filter(diff_men_women>0.1)
tab2<-data.frame(men=t(d1 %>% summarise_at(vars(men_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 women=t(d1 %>% summarise_at(vars(women_logit),funs(eff,se,p,sup20,sup80,fact,n))))
d1 <- d %>% filter(diff_inc90_inc10>0.08 & topic=="Economy")
tab3<-data.frame(men=t(d1 %>% summarise_at(vars(men_logit),funs(eff,se,p,sup20,sup80,fact,n))),
                 women=t(d1 %>% summarise_at(vars(women_logit),funs(eff,se,p,sup20,sup80,fact,n))))
tab <- cbind(tab1,tab2,tab3)
tab1 <- data.frame(t(tab))
row.names(tab1) <- c("Men",
                     "Women",
                     "Men ",
                     "Women ",
                     "Men  ",
                     "Women  ")
kable(tab1[-3],booktabs=T, align="c",linesep = '\\addlinespace', col.names = c("Effect (logit coefficient)","Standard error", "Predicted probability of policy change at 20% support","Predicted probability of policy change at 80% support","Relative change in probability","N"),
      caption = "Policy responsiveness by gender") %>%
  column_spec(1,width="1.4in") %>%
  column_spec(2,width="0.5in") %>% 
  column_spec(3,width="0.55in") %>% 
  column_spec(4:5,width="0.8in") %>% 
  column_spec(6,width="0.5in") %>% 
  column_spec(7,width="0.4in") %>% 
  pack_rows("All issues",1,2) %>%
  pack_rows("When preferences diverge",3,4) %>%
  pack_rows("Economic policy, preferences diverge between rich and poor",5,6) %>%
  footnote("Bivariate logistic regression models (rows). The dependent variable is a dichotomous measure of whether or not the policy change was adopted within four years of the time of the survey question. Rows 3 and 4 based on the set of issues where the preference distance between men and women is larger than 10 points. Rows 5 and 6 based on the set of issues where the preference distance between the 10th and 90th income percentile is larger than 8 points (cf. Table A3).*p<0.1; **p<0.05; ***p<0.01 *p<0.1; **p<0.05; ***p<0.01",threeparttable = T) %>%
  kable_styling(latex_options = "HOLD_position")


# Table H1 
## Access to Norwegian Election Studies (NES) can be ordered for from https://www.nsd.no/nsddata/serier/norske_valgundersokelser.html
## NES 1965
Valg1965 <- read_dta("NES 1965.dta")
Valg1965$women <- car::recode(Valg1965$v264, "2=1; 1=0; else=NA")
Valg1965$foreign <- car::recode(Valg1965$v294,"8=1; 1:4=0; else=NA")
Valg1965$low_edu <- car::recode(Valg1965$v286,"1:2=1; 3:9=0; else=NA")
Valg1965$married <- car::recode(Valg1965$v266,"1=1; 2:4=0; else=NA")
Valg1965$age <- as.numeric(Valg1965$v265)
## NES 1989
Valg1989 <- read_dta("NES 1989.dta")
Valg1989$women <- car::recode(Valg1989$V16, "2=1; 1=0; else=NA")
Valg1989$foreign <- car::recode(Valg1989$V72,"88=1; 1:20=0; else=NA")
Valg1989$low_edu <- car::recode(Valg1989$V528,"1:4=1; 5:9=0; else=NA")
Valg1989$married <- car::recode(Valg1989$V468,"1=1; 2:5=0; else=NA")
Valg1989$age <- 89-Valg1989$V14
## NES 2013
Valg2013 <- read_dta("NES 2013.dta")
Valg2013$women <- car::recode(Valg2013$Kjonn, "2=1; 1=0; else=NA")
Valg2013$foreign <- car::recode(Valg2013$Bak30a,"2=1; 1=0; 8:9=0; else=NA")
Valg2013$low_edu <- car::recode(Valg2013$Bak28,"0:21=1; 22:99=0; else=NA")
Valg2013$married <- car::recode(Valg2013$Bak1,"1=1; 2:9=0; else=NA")
Valg2013$age <- Valg2013$Alder
## Make data table
prop <- function(x)round(prop.table(table(x))[2]*100,1)
demo <- data.frame(Statistic=c("% Women",
                               "% With only primary, lower secondary, or no education",
                               "% Married",
                               "% Born/raised abroad",
                               "Mean age"),
                   `1965`=c(prop(Valg1965$women),
                            prop(Valg1965$low_edu),
                            prop(Valg1965$married),
                            prop(Valg1965$foreign),
                            round(mean(Valg1965$age,na.rm=T),1)),
                   `1989`=c(prop(Valg1989$women),
                            prop(Valg1989$low_edu),
                            prop(Valg1989$married),
                            prop(Valg1989$foreign),
                            round(mean(Valg1989$age,na.rm=T),1)),
                   `2013`=c(prop(Valg2013$women),
                            prop(Valg2013$low_edu),
                            prop(Valg2013$married),
                            prop(Valg2013$foreign),
                            round(mean(Valg2013$age,na.rm=T),1)) )
## Make presentable table 
kable(demo,booktabs=T, align="l",linesep = '\\addlinespace',
      caption = "Survey demography",
      col.names = c("Statistic","1965","1989","2013")) %>%
  column_spec(1,width="2in") %>%
  footnote("Source: Norwegian National Election Studies.",threeparttable = T) 


# Table H2 and H3: See seperate file called "Code for Appendix Tables H2 and H3".


###################### END #############################




